import matplotlib.pyplot as plt
import numpy as np
import glob


plt.rc('axes', titlesize=14)     # fontsize of the axes title
plt.rc('axes', labelsize=14)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=14)    # fontsize of the tick labels
plt.rc('ytick', labelsize=14)    # fontsize of the tick labels

fl=glob.glob('currentGP-TSD/**/*.osc')
print(fl)
#_3to1left=np.loadtxt('3to1left.dat')
#_3to1right=np.loadtxt('3to1right.dat')
#slope3to1, intercept3to1 = np.polyfit(_3to1left[:,0],_3to1left[:,1],1)

iii=[]
for f in fl:
    indx, a,e,i=np.loadtxt(f,unpack=True,usecols=(0,2,3,4))
    ind=np.argwhere( (a>2.45) & (e<0.15) & (i<15))
    print(f, indx[ind], a[ind], e[ind], i[ind])
    
    ind=np.argwhere(i<15)
    plt.scatter(a[ind],e[ind],c='#888888',s=5)

#exit(0)

Athor=[2.38,  0.1,  0.15]
#AthorReg=(Athor[0]-0.1, Athor[0]+0.1)
AthorReg=(2.3, 2.48)

x=np.linspace(2,2.5,100)
y=1-1.8/x
plt.plot(x,y,c='k')
#y=1-2.0/x
#plt.plot(x,y,c='#BBBBBB')


x=np.linspace(AthorReg[0],AthorReg[1])
yUp=1-2.0/x
yDo=np.zeros(len(x))
plt.fill_between(x, yUp, yDo, color='#BBBBBB', alpha=0.5)

plt.scatter(Athor[0],Athor[1],s=30,c='k',zorder=10)

#x=np.linspace(2,2.5,100)
#plt.fill_between(x, np.ones(len(x)), -12.5*(x-2.5), color='#BBBBBB', alpha=0.5)
plt.ylim(0,1.0)
plt.xlim(2.0, 2.5)
#plt.xlabel('Semimajor axis, $a$ (au)')
plt.ylabel('Eccentricity, $e$')
plt.text(2.364, 0.055, "Athor")

ax = plt.gca()
ax.axes.xaxis.set_ticklabels([])

plt.tight_layout()
plt.savefig('Fig2A.pdf',dpi=300)
